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ABSTRACT 

Tseliakhovich & Hirata recently discovered that higher-order corrections to the cosmological 
linear-perturbation theory lead to supersonic coherent baryonic flows just after recombination 
(i.e. z ps 1020), with rms velocities of ~30 km/s relative to the underlying dark-matter distri- 
bution, on comoving scales of < 3 Mpc h~ l . To study the impact of these coherent flows we 
performed high-resolution N-body plus SPH simulations in boxes of 5.0 and 0.7 Mpc ft. -1 , 
for bulk-flow velocities of (as reference), 30 and 60 km/s. The simulations follow the evo- 
lution of cosmic structures by taking into account detailed, primordial, non-equilibrium gas 
chemistry (i.e. H, He, H2, HD, HeH, etc.), cooling, star formation, and feedback effects from 
stellar evolution. We find that these bulk flows suppress star formation in low-mass haloes 
(i.e. M vir < 1O 8 M until z ~ 13), lower the abundance of the first objects by — 1% - 20%, 
and as consequence delay cosmic star formation history by ~ 2 x 10 7 yr. The gas fractions 
in individual objects can change up to a factor of two at very early times. Coherent bulk flow 
therefore has implications for (i) the star-formation in the lowest-mass haloes (e.g. dSphs), 
(ii) the start of reionization by suppressing it in some patches of the Universe, and (iii) the 
heating (i.e. spin temperature) of neutral hydrogen. We speculate that the patchy nature of 
reionization and heating on several Mpc scales could lead to enhanced differences in the HI 
spin-temperature, giving rise to stronger variations in the HI brightness temperatures during 
the late dark ages. 
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1 INTRODUCTION 

Cosmic stru cture-formation models are based on Jeans theory 
jjeanslll902h applied t o primordial matter fl uctuations in an ex- 
panding Universe (e.g. iKomatsu et al]|2010f) . The matter density 
perturbations are linearly expanded to first order, which is supposed 
to work w ell as long as the density contrast is much smaller than 
unity (e.g. |Peebleslll974l : IColes & Lucchinll2002l ; ICiardi & Ferraral 
2005), and the equations that govern the evolution of the dark and 
baryonic matter become linear. To study highly non-linear evolu- 
tion when the density contrast approaches and far exceeds unity, 
however, fully numerical simulations are required. Accordin g to 
the concordance cosmological model (e.g. K omatsu et alj2010h the 
Universe is flat and has an expansion rate of Hq ~ 71 km/s/Mpc. 
The cosmological constant (or dark energy) represents the main 
energy-density content, £Io,a — 0.7, whereas dark matter con- 
tributes f2o,M — 0.26 and baryonic matter is only a small fraction 
of f2o,b — 0.04. Within this context, baryonic structures arise from 
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in-fall and condensation of gas into growing dark-matter potentia l 
wells (e.g. lGunn & Gottll972| ; |Peebleslll974l ; IWhite & Reeslll978h . 
Both in-fall and condensation of gas only take place after the time 
of recombination (z ~ 1020), when baryons and photons decouple 
dKomatsu et alj2oToh . On the other hand, dark matter perturbations 
can already start growing when their densit y exceeds that of radia- 
tion (i .e. at z cq ~ 3100 — 3200). Recently Tseliakhovich & Hiratal 
d2010h for the first time included quadratic terms in the cosmo- 
logical perturbation theory to account for the advection of small- 
scale perturbations by large-scale velocity flows. They find that, at 
the redshift of decoupling z ~ 1020, coherent supersonic flows 
of the baryons relative to the underlying dark-matter distribution 
are formed on scales of a few Mpc or less, with typical veloci- 
ties on those scales with a rms value of ~ 30 km/s, sourced by 
density perturbations up to scales of ~ 100 Mpc. These gas bulk 
velocities could actually suppress the formation of the very early 
and low-mass stru ctures, induce higher baryon acoustic oscillations 
(Dalai et al. 2010), and have further implications on high-redshift 
galaxy clustering, 21-cm studies, and reionization. In fact, as we 
find, they might impact the formation of the lowest-mass galax- 
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ies (e.g. dSph) on scales as large as the local group, depending 
of the magnitude of the bulk flow at decoupling. The effect can 
thus cause a spatially varying bias in baryonic structure forma- 
tion, whereas small-scale baryonic structures are suppressed more 
in regions where the bulk-velocity is larger. Detailed numerical 
simulations considering this non-linear effect are not yet available 
and the very high dynamic range required is not foreseen in the 
near future. It therefore looks very challenging to quantitatively 
assess these effects in the hi ghly non-linear regime discussed by 
Tseliakhovich & Hirata (2010). Indeed, one would need to resolve, 
at the same time, large scales (of the order of the horizon) and ex- 
tremely small scales (of the order of the Jeans length). Given these 
extreme difficulties, we adopt the following approach: we focus on 
small scales where the effect is expected to be largest, at least prop- 
erly dealing with all the detailed physics and chemistry, and we 
assume a co nstant bulk velocity shift in the initial conditions as de- 
termined bv lTseliakhovich & Hiratal d2010h on the relevant scales. 
Even though this approximation is still crude, we believe it is a rea- 
sonable assumption to study the lowest-mass structures in the box 
(which have virial radii much smaller than the box size), we expect 
it to be able to provide the first quantitative results in the highly 
non-linear regime and assess its effects on early structure forma- 
tion. In this paper, we make a first step towards understanding the 
effects of coherent supersonic bulk flow on the formation of the 
smallest baryonic structures, in particular focusing on the suppres- 
sion of star formation and the delay of reionization. In Sect. [2] we 
describe our simulations, and in Sect. [3] we show our results. We 
summarize and conclude in Sect. [4] 



2 SIMULATIONS 

We perform nu merical simulati ons by using the parallel tree/SPH 
P-Gadget2 code (Springel 2005), which implements gravity, hydro- 
dynamics, molecular and atomic evolution, cooling, population III 
and population II star formation, supernova/wind feedback, and a 
full chemistry network involving e _ , H, H + , H , He, He + , He ++ , 
H 2 , Hj~, D D + , HP, and HeH + (for more deta ils see lMaio et~al1 
l200dl200l[20ld : lMaidl2009l ; lMaio et alj|20ld) . Particularly im- 
portant for our purposes is the capability of the code of following 
gas collapse down to the catastrophic molecular cooling branch, by 
resolv ing the Jeans sca les of primordial structures with ~ 10 2 par- 
ticles (Maio et al.ll2009h . and completely ful filling the correct reso- 
lution requireme nts for SPH simulations (see Ba te & Burkertl 19971 ; 
Mai o et al J2009L for details). We generate the initial conditions by 
using the N-GenIC code and add bulk velocity shifts for the gas in 
the x-direction i>b, x = 0, 30, 60 km/s (referred to as HR-ShiftOO, 
HR-Shift30, and HR-Shift60, respectively) at recombination epoch 
and p roperly scaled down (according to iTseliakhovich & Hiratal 
2010) to the simulation initial redshift of Zi n = 100. The value 
Wb,x ~ 60 km/s has been considered to account for statistical vari- 
ations from the expected rms value of ~ 30 km/s, but it should 
be considered as an upper limit which can be found in less than 
1% of the volume of the Universe. Even though we expect that the 
presence of bulk flows at higher redshift does not affect the struc- 
ture formation process, we nevertheless run additional simulations 
with z in = 1020, with v b , x = 0,60 km/s (HR-ShiftOORec, HR- 
Shift60Rec). We should caution the reader that the Gadget code 
implements N-body/SPH structure formation, cooling, star forma- 
tion, and feedback effects in matter-dominated universes, while 
contributions from radiation or from the primordial hot plasma 
in the cosmological dynamics are neglected. Therefore, its use at 



z » 200 might be questionable. In any case, the results are very 
similar to the simulations with Zi n = 100, and thus we don't dis- 
cuss them any further. Our reference simulations have a box side 
of 0.7Mpc//i ~ 1 Mpc, and are run in the frame of the stan- 
dard ACDM cosmological model (fio.A = 0.7, fio.M = 0.26, 
fi ,b = 0.04, h = 0.7, a 8 = 0.9, n = 1). Matter and gas fields are 
sampled with 320 3 particles, respectively, allowing a resolution of 
~ 1O 2 M /h for the gas component and ~ 7.5 x 1O 2 M /h for the 
d ark-matter component. T he chemical abundances are initialized as 
in lMaio et alj d2007Ll2010b . To check if the effects of primordial su- 
personic gas flows on structure formation could be seen in bigger 
but lower-resolution simulations, we also run a set of three larger 
boxes (LR-ShiftOO, LR-Shift30, and LR-Shift60) with a size of 
5 Mpc//j, Wb, x = 0, 30, 60, km/s, and with the same parameters de- 
scribed before. This box is comparable to the smallest scale where 
the bu lk flows are expected according to Tseliakhovich & Hiratal 
(2010). The corresponding mass resolution is ~ 4x lO 4 M //iand 
~3x 10 5 Mq /h, for gas and dark-matter species. In this case, we 
did not find any relevant differences among the three runs, implying 
that the effect is observable only if the numerical resolution is high 
enough to resolve halo masses of ~ lO 4 M /h. This highlights the 
difficulty of current large-scale simulations to address these effects. 
In TableQ]we list the main features of the numerical set-ups. 



3 RESULTS 

We start by looking at the overall evolution of the star formation 
rate density of the Universe, and then we consider the main statis- 
tical quantities involved in the cosmological evolution process, i.e. 
the distributions of dark-matter haloes and gas clouds. 



3.1 Star formation rate density 

The behavior of the star formation rate densities is plotted in Fig. 
[T] for all the cases considered before. What emerges is a system- 
atic shift in the onset of star formation, with a delay increasing for 
higher values of «b, x . For the «b, x = Okm/s case, the onset hap- 
pens at z ~ 16.3, when the Universe is roughly 2.5 x 10 8 yr old, 
while in the Ub, x = 30 km/s and 60 km/s cases, star formation 
sets in at z ~ 16.1 and 15.6, respectively, corresponding to a delay 
of ~ 10 7 yr and ~ 2 x 10 7 yr compared to the case without ve- 
locity shift. Differences at very early times are almost one order of 
magnitude and at later times a factor of a few, persisting down to 
z ~ 13, when the global trends rapidly converge to similar values. 
The reason why this happens is due to the higher kinetic energy 
given to the gas at the recombination epoch: because of that, first 
dark-matter haloes can retain less material, gas condenses more 
slowly hindering molecule formation, the corresponding cooling is 
less efficient, and the resulting star formation is delayed. In partic- 
ular, this is very relevant for small pri mordial haloes, with masse s 
of ~ 1O 4 M - 10 8 M o (see also ITseliakhovich & Hirat3l201Ch . 
whose dimensions are comparable to or smaller than the baryon 
Jeans length. In this case in fact, the gas cannot fragment within 
the halo and partially or entirely flows out of it, since larger dark- 
matter potential wells are needed to retain the gas. Therefore, this 
effect on small scales influences the growth of larger objects, as 
well. From this we conclude that the inflow and condensation of 
gas, and thus also star formation, is already strongly suppressed in 
haloes in the mass-range < 10 s M before reionization. Those low- 
mass halos surviving to the present day, or even those that merge 
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Table 1. Parameters adopted for the simulations. 
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Figure 1. Total star formation rates for bulk velocity shifts of km / s (black 
solid line), 30 km/s along the x-axis (blue dotted line), 60 km/s along the 
x-axis (red dashed line) are shown for the simulation box of side ~ 1 Mpc. 

into larger masses, would therefore have a lower baryonic and stel- 
lar mass fraction. We further note that the effect of suppressing star 
formation is similar to that of reionization, but it has a very different 
physical origin and occurs well before reionization itself. 

3.2 Dark-matter haloes and gas clouds 

In each simulation, we identify cosmic structures with the use of 
a friend-of-friend (FoF) algorithm. Each object is represented by 
all those particles (dark-matter, gas and star) closer than a limit- 
ing linking length of 20% the mean inter-particle separation. The 
dark-matter halo and gas cloud distributions are shown in Fig. [2] 
at different redshift as indicated by the legends, for the runs start- 
ing at Zi n — 100. In general, at very high redshift we cannot give 
solid conclusions, since small-number statistic dominates the sam- 
ple. At later times (i.e. z < 23), some systematic behavior clearly 
emerges. The mass distributions of dark matter haloes range be- 
tween ~ 1O 4 M0//i and lO 7 M //i. The three cases of «b,x = 
0, 30, 60 km/s show a slight decrement with increasing Ub, x , but 
only at a few percent level. Hence, as expected, the effect of the 
bulk-flow of gas on the overall growth of the dark-matter dominated 
haloes is minimal. More clear differences arise from the analyses 
of the gas clouds. In Fig[2] we also show the mass distribution of 



the gaseous component of the objects and we clearly see decre- 
ments up to tens of percent, mostly evident in the «b, x = 60 km/s 
case. Interestingly, there is a cascade effect from smaller to larger 
masses (expected from the simple analyses of the star formation 
rates), which delays and slightly hinders the accumulation of gas in 
the dark-matter haloes (i.e. suppression of gas inflow in lower-mass 
haloes also suppresses the gas accumulation in high-mass haloes 
through accretion of low-mass haloes). As expected, the average 
gas cloud mass decreases for larger t>b,x- At z ~ 23, it drops from 
~ 1.70xl0 3 M Q //i(u b , x = km/s) down to ~ 1.62 x lO 3 M //i 
Ob,x = 30 km/s), and ~ 1.39 x lO 3 M //i (« b , x = 60 km/s). 
While at z ~ 15, it is ~ 4.41 x lO 3 M //i (t> b , x = Okm/s), 
~ 4.26 x 10 3 M Q //i Cu b , x = 30 km/s), and ~ 4.06 x W 3 M e /h 
(fb,x = 60 km/s). The differences correspond to decrements with 
respect to the reference run of ~ 6% — 3% (fb,x = 30 km/s) and 
~ 20% - 10% («b, x = 60 km/s) at z ~ 23 - 15, respectively. 
Also the total abundance of gas clouds shows a similar trend, with 
mass-weighted number counts dropping, with respect to the refer- 
ence run, by ~ 5% - 3% (v b ^ = 30 km/s) and ~ 20% - 50% 
(Wb, x = 60 km/s) at z ~ 23 — 15. We have verified that the ex- 
cess small clumps observed in the low-mass tail of the distribu- 
tions (at masses < 5OOM //i) for «b, x > are due to residual 
gas which, differently from the reference run, is not trapped by 
the dark-matter haloes and, later collapses to form such clumps. 
It should be noted though that statistical errors in this mass regime 
are substantial, due to the small number of particles contained in 
a halo. As a more specific example, in Fig. [3] we show the evolu- 
tion of the gas content (top panel) and gas fraction (bottom panel) 
of the first (and most massive) halo formed in the simulations. The 
gas mass ranges between ~ 1O 3 M //i and lO 7 M //i, and the 
corresponding dark-matter mass is about one order of magnitude 
larger. The gas mass variations due to the different values of Db,x 
are significant at z > 16, while at smaller redshift the total amount 
converges to ~ 2 — 6 x 10 M jh. In this respect, the gas fractions 
in the halo are more informative, and suggest a depletion of gas for 
higher «b, x up to a factor of 2. Indeed, at z > 24, the gas fraction 
for the Ub. x = km/s case is ~ 0.04, while in the Ub. x = 30 km/s 
and Ub. x = 60 km/s cases it drops to ~ 0.03 and below 0.02, re- 
spectively. Similar strong variations are observed also at later times 
(e.g. at z ~ 23), but the general trend converges at z < 16, when 
the discrepancies among the cases become < 10% and all the gas 
fractions reach values of ~ 0.10. The same effects are observed in 
the gas mass and density profiles as functi on of the (physical) ra- 
dius - see Fig.|4]and compare to Fig. 2 in lTseliakhovich & Hiratal 
At early times (left panels) the amount of gas falling in the 
halo depends on «b, x . and the larger u b ,xi the more gas is depleted 
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Figure 2. Top row: Dark-matter halo distributions at different redshifts (see labels) for the 1 Mpc side boxes. Bottom row: Gas cloud distributions at different 
redshifts (see labels) for the 1 Mpc side boxes. Histograms refer to velocity shifts of km/s (black crosses), 30 km/s along the x-axis (blue rhombi), 60 km/s 
along the x-axis (red triangles). 
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Figure 3. Evolution of the gas content (top) and of the gas fraction (bottom) 
for the first (and most massive) halo formed in the simulations. Data refer 
to velocity shifts of Okm/s (solid lines, black crosses), 30 km/s (dotted 
lines, blue rhombi), 60 km/s (dashed lines, red triangles). 

in the distribution tails. At larger radii, the mass depletion is up to 
~ 50% and smoothly decreases towards the core, where the gas is 
even denser and more abundant. This can be seen as a consequence 
of the fact that there is less gas falling in, and it sinks more easily 
in the potential wells of the halo. Of course, given the low number 
of SPH particles reaching those regimes, the statistical errors might 
be quite significant and ad hoc simulations are needed to address 
this issue in more detail. 



4 SUMMARY AND CONCLUSIONS 

Stimu lated by very recent analytical work jTseliakhovich & Hiratal 
2010), we have run cosmological, high-resolution numerical N- 
body/SPH simulations including detailed, non-equilibrium chem- 
istry evdutic^coohng^^ar formation and feedback effects (see 
e.g. lMaio et alJ[2007L l2010h . which consistently follow the gas be- 
havior down to the catastrophic cooling regime, and correctly re- 
solve the small, primordial Jeans scales with ~ 10 2 SPH parti- 
cles. To take into account very early, supersonic, coherent, Mpc- 
scale flows, generated at recombination time from the advection of 
small-scale perturbations by large-scale velocity flows, we set ini- 
tial gas velocity shifts along the 2>direction, u b x = 0, 30, 60 km/s 
based on analytic estimates of the rms differences in dark-matter 
and baryonic gas velocities at z = 1020. We find that in the simula- 
tions with larger t>b,x the onset of star formation is delayed by some 
~ 10 7 yr, and it is smaller by a factor of a few in comparison to the 
Wb,x = 0km/s case, until redshift z ~ 13. As a consequence, 
the beginning of the reionization process is delayed by the same 
amount and initially driven by more massive objects. More ad hoc 
calculations should be done to quantify the differences in terms of 
evolution, topology and photon budget. This might also have con- 
sequences in terms of detection of the 21cm signal from neutral 
hydrogen, because th e building up of a Lya background (see e.g. 
iRobertson et al.l2010l for a recent review) ne cessary for the observ- 
ability of the line (e.g. lFurlanetto et al]|2006l) would be delayed as 
well. We speculate that this could lead to more patchy reionization 
and HI heating, with larger HI brightness-temperature fluctuations. 
More research is needed to quantify these effects. Because we find 
that the amount of gas in the first haloes is expected to drop signif- 
icantly (up to a factor of 2 at early times, i.e. ~ 50%), there could 
also be interesting consequences in terms of feedback effects, as 
e.g. less gas in small mass objects should make it easier to fur- 
ther remove gas via photoevaporation, winds or SN explosions. It 
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Figure 4. Gas profiles as a function of the (physical) radius at redshift z = 23 (left panels) and z = 15 (right panels) for the data referring to the simulations 
with velocity shifts of km/s (solid black lines), 30 km/s (dotted blue lines) and 60 km/s (dashed red lines). Top panels show the gas mass profiles, while 
bottom panels show the gas density profiles. 



will be interesting to investigate in more details the interplay be- 
tween bulk-flow induced feedback effects. Slight differences (at 
~ 1% level) in the dark-matter statistical distributions also arise, 
but stronger differences (up above ~ 10% level) are seen in the 
gas clouds. A cascade effect transfers the inability of primordial, 
small haloes to completely retain the gas to larger, lower-redshift 
haloes, whose gaseous components are systematically smaller for 
«b,x > km/s. In this way, the original hindering for the very first 
objects is transferred to the structure which form later on, affect- 
ing their growth history. Numerical studies performed on lower- 
resolution simulations show no detectable differences among the 
different cases, and this means that the lack of numerical resolution 
can be a strong limitation to probe higher-order corrections to lin- 
ear perturbation theory, and SPH particle masses of ~ 10 2 M o //i 
are required. The precise initial redshift does not alter significantly 
our findings either. 

Our general conclusions from this albeit preliminary, but first fully 
non-linear, study of the effects of bulk gas flows on non-linear cos- 
mic structure growth are that it: (i) delays early star formation, 
which affects both reionization (by delaying it to lower redshifts) 
and the heating of gas of higher redshift, possible affecting to- 
tal/global emission and absorption features of HI against the CMB 
by making these effects more patchy; (ii) suppresses star formation 
in the lowest-mass haloes over all redshifts which has a cascading 
effect over nearly 300 million years on higher-mass haloes which 
are their merger-products. Because these small early objects, if they 
survive till the present day, are expected to be the parent population 
of e.g. dSph galaxies, the effct of the bulk flows can also have an 
impact on the stellar/gas content of the lowest-mass dwarf satel- 
lites around more massive galaxies (e.g. the MW) and be highly 
spatially dependent (i.e. the effect is only present where there were 
large bulk flows at decoupling). Th is could result in a strongly spa- 
tiall y varying bias as s uggested by Tseliak hovich & Hiratal d201Ch 
and Dalai et a I] d2010h . In addition, suppression of gas infall and 
condensation on these small scales coul d play an impor tant role in 
the "missing satellite problem" (see e.g. lKravtsowOlOl and refer- 
ences therein). 
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